//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : InteractiveWaves.cpp
//----------------------------------------------------------------------------

#include <iostream>
#include <algorithm> // min + max
#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "stdio.h"
#include "InteractiveWaves.h"
#include "FloatUtils.h"
#include "StringUtils.h"
#include "FileReader.h"
#include "WavesComputer.h"
#include "Solver.h"
#include "assertplus.h"

#define MAXIMUM_DEPTH 50

#ifdef _WINDOWS
#define _CRT_SECURE_NO_WARNINGS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#pragma warning(disable: 4305) // truncation from 'double' to 'float'
#pragma warning(disable: 4056) // overflow in floating-point constant arithmetic
#pragma warning(disable: 4756) // overflow in constant arithmetic
#pragma warning(disable: 4996) // This function or variable may be unsafe
#endif

static double getBilinearInterpolation(const Array2D<REAL>& a_waves,     ///< [m] wave heights (with solver ghosts)
                                       double a_dx,                      ///< [m] interval in direction 1
                                       double a_dy,                      ///< [m] interval in direction 2
                                       const mml::Vector2d& a_position); ///< [m] x,y position of point of interest

/// constructor
InteractiveWaves::InteractiveWaves()
: m_pWC(NULL)
, m_ph0(0)
, m_pU1(0)
, m_pU2(0)
, m_pwet(0)
, m_pzeta(0)
, m_pzetam1(0)
, m_pphi(0)
, m_pphim1(0)
, m_pzetaNet(0)
, m_Nx1(0)
, m_Nx2(0)
, m_dx1(0)
, m_dx2(0)
, m_g(0)
, m_h0_c(0)
, m_U1_c(0)
, m_U2_c(0)
, m_wvin_spec_count(0)
, m_wvin_spec_ampl()
, m_wvin_spec_phase()
, m_wvin_spec_freq()
, m_wvin_spec_dir()
, m_steps(0)
, m_skipani(0)
, m_h0_min(1e300)
, m_h0_max(0)
, m_save(false)
{
}

/// destructor
InteractiveWaves::~InteractiveWaves()
{
    delete m_pWC;
    delete m_ph0;
    delete m_pU1;
    delete m_pU2;
    delete m_pwet;
    delete m_pzeta;
    delete m_pzetam1;
    delete m_pphi;
    delete m_pphim1;
    delete m_pzetaNet;
}


/// reads parameters from datafile(s) and creates WavesComputer
bool InteractiveWaves::setup(const char *parameterfile, ///< file with level, waves, currents
                             const char *solvername, ///< name of type of solver to use
                             REAL dt, ///< [s] delta T
                             REAL pstart) ///< [s] fade-in time of pulses and incoming waves
{
    m_g = getGravity();
    VariationalModel varmod = undefVM;
    REAL kappa_c = 0, tau_brk = 0;
    bool short_wave = true, upwind = false;
    REAL bcN = 0, bcE = 0, bcS = 0, bcW = 0, somfeld = 0, k_wrbc = 0;
    REAL eps = 0, lp = 0, psitol = 0;
    Initialization init = noInit;
    REAL time = 0.0;
    std::string file = "";
    try
    {
        FileReader fr(parameterfile);
        do
        {
            std::string section = fr.ReadString();
            ASSERT(section != "[Physics]", "verwijder Physics sectie uit par-file"); // TODO remove
            if (section == "[WaveModel]")
            {
                varmod = VariationalModel(fr.ReadInt());
                kappa_c = fr.ReadFloat();
                tau_brk = fr.ReadFloat();
                short_wave = fr.ReadBoolean();
                upwind = fr.ReadBoolean();
            }
            else if (section == "[Area]")
            {
                setupArea(fr, dt, &init, &time, &file);
            }
            else if (section == "[Boundaries]")
            {
                bcN = fr.ReadFloat();
                bcE = fr.ReadFloat();
                bcS = fr.ReadFloat();
                bcW = fr.ReadFloat();
                somfeld = fr.ReadFloat();
                k_wrbc = fr.ReadFloat();
            }
            else if (section == "[IncomingWaves]")
            {
                setupIncomingWaves(fr);
            }
            else if (section == "[Computations]")
            {
                eps = fr.ReadFloat();
                lp = fr.ReadFloat();
                psitol = fr.ReadFloat();
            }
        }
        while (!fr.Eof());
        m_pWC = new WavesComputer(m_Nx1, m_Nx2, *m_ph0, *m_pU1, *m_pU2, *m_pwet,
                                  m_pzeta, m_pzetam1, m_pphi, m_pphim1, init, time, dt,
                                  pstart, solvername, m_dx1, m_dx2, m_g,
                                  tau_brk, m_wvin_spec_count, m_wvin_spec_ampl,
                                  m_wvin_spec_phase, m_wvin_spec_freq, m_wvin_spec_dir,
                                  int(bcN), int(bcE), int(bcS), int(bcW), int(somfeld), k_wrbc, varmod,
                                  kappa_c, psitol, short_wave, upwind, eps, int(lp));
        if (init == fileInit)
        {
            m_pWC->load(file.c_str());
        }
        return true;
    }
    catch (...)
    {
        return false;
    }
}


/// computes waves in one simulation step
void InteractiveWaves::compute(REAL dt, ///< delta T
                               Array2D<REAL> &pulseField, ///< array with pulse pattern
                               Array2D<REAL> **ppHeight, ///< resulting heightfield
                               Array2D<REAL> **ppNetHeight, ///< resulting heightfield wihout pulse 'hole'
                               Array2D<REAL> **ppSpeed, ///< resulting speed potential
                               double *pCgTime) ///< time used by solver, for testing purposes
{
    // compute m_pzeta and m_pphi
    m_pWC->compute(dt, pulseField, pCgTime);

    // compute m_pzetanet
    correctForPulses(pulseField);

    // save wave field if desired
    if (m_save)
    {
        m_pWC->save("zetaphi.wave");
        m_save = false;
    }
    
    // return values
    *ppHeight = m_pzeta;
    *ppNetHeight = m_pzetaNet;
    *ppSpeed = m_pphi;
}


/// returns parameters of simulated area
void InteractiveWaves::getArea(int *pNx1, ///< number of points in x1 direction
                               int *pNx2, ///< number of points in x2 direction
                               REAL *pdx1, ///< grid size in x1 direction (m)
                               REAL *pdx2, ///< grid size in x2 direction (m)
                               REAL *pg) ///< gravity
const
{
    ASSERTPOINTER(pNx1);
    ASSERTPOINTER(pNx2);
    ASSERTPOINTER(pdx1);
    ASSERTPOINTER(pdx2);
    ASSERTPOINTER(pg);
    *pNx1 = m_Nx1;
    *pNx2 = m_Nx2;
    *pdx1 = m_dx1;
    *pdx2 = m_dx2;
    *pg = m_g;
}


/// returns bathymetry data
void InteractiveWaves::getBathymetry(Array2D<REAL> **ppField, ///< heightfield
                                     REAL *pMin, ///< minimum height
                                     REAL *pMax, ///< maximum height
                                     Array2D<SHORT> **ppWet) ///< wet/dry field
const
{
    ASSERTPOINTER(ppField);
    ASSERTPOINTER(pMin);
    ASSERTPOINTER(pMax);
    ASSERTPOINTER(ppWet);
    *ppField = m_ph0;
    *pMin = m_h0_min;
    *pMax = m_h0_max;
    *ppWet = m_pwet;
}


/// returns simulator time
REAL InteractiveWaves::getTime() const
{
    return m_pWC->getTime();
}

/// saves wave height and speed potential to file
void InteractiveWaves::save()
{
    m_save = true;
}

//----------------------------------------------------------------------------
// function InteractiveWaves::getHeight
//----------------------------------------------------------------------------
/// determines wave height at specified location
/** @return wave height (z value) */
//----------------------------------------------------------------------------
REAL InteractiveWaves::getHeight(mml::Vector2d const& a_xy) const
{
	return getBilinearInterpolation(*m_pzeta, m_dx1, m_dx2, a_xy);
}


//----------------------------------------------------------------------------
// function InteractiveWaves::getNetHeight
//----------------------------------------------------------------------------
/// determines net wave height at specified location
/** @return net wave height (z value) */
//----------------------------------------------------------------------------
REAL InteractiveWaves::getNetHeight(mml::Vector2d const& a_xy) const
{
	return getBilinearInterpolation(*m_pzetaNet, m_dx1, m_dx2, a_xy);
}


//----------------------------------------------------------------------------
// function InteractiveWaves::getSpeed
//----------------------------------------------------------------------------
/// determines speed potential at specified location
/** @return speed potential */
//----------------------------------------------------------------------------
REAL InteractiveWaves::getSpeed(mml::Vector2d const& a_xy) const
{
	return getBilinearInterpolation(*m_pphi, m_dx1, m_dx2, a_xy);
}

#define BUFLEN 200
static void ReadLine(FILE *fp, REAL *a, REAL *b, REAL *c, bool *d)
{
    char buffer[BUFLEN];
    char *words[4];
    char *ptr;
    char* ignore;
    ignore = fgets(buffer, BUFLEN, fp);
    ptr = StringUtils::trim(buffer);
    StringUtils::split(ptr, words, 2, '%');
    StringUtils::split(words[0], words, 4, ' ');
    *a = atof(StringUtils::trim(words[0]));
    *b = atof(StringUtils::trim(words[1]));
    *c = atof(StringUtils::trim(words[2]));
    *d = strcmp(StringUtils::trim(words[3]), "T") == 0;
}

//----------------------------------------------------------------------------
// function getBilinearInterpolation
//----------------------------------------------------------------------------
/// interpolates wave height at specified x,y coordinate
/** @return z coordinate of wave surface */
//----------------------------------------------------------------------------
static
double getBilinearInterpolation(const Array2D<REAL>& a_waves_1, ///< [m] wave heights (with solver ghosts)
                                double a_dx,                     ///< [m] interval in direction 1
                                double a_dy,                     ///< [m] interval in direction 2
                                const mml::Vector2d& a_position) ///< [m] x,y position of point of interest
{
    // determine indices of centroid location in waveheight array
    double x, y;
    modf(a_position(0) / a_dx, &x);
    modf(a_position(1) / a_dy, &y);
    int i = int(x);
    int j = int(y);
    // cannot interpolate outside array!
	if (i < 0 || i >= int(a_waves_1.dim1()) - 3 || j < 0 || j >= int(a_waves_1.dim2()) - 3)
	{
	    return 0.0;
	}

    // interpolate waveheight at position
    // see http://en.wikipedia.org/wiki/Bilinear_interpolation
    x = a_position(0);
    y = a_position(1);
    REAL Q11 = a_waves_1[i     + 1][j     + 1]; // waves array has ..
    REAL Q12 = a_waves_1[i     + 1][j + 1 + 1]; // .. border of ghost points ..
    REAL Q21 = a_waves_1[i + 1 + 1][j     + 1]; // .. for solver ..
    REAL Q22 = a_waves_1[i + 1 + 1][j + 1 + 1]; // .. so add 1 to indices
    double x1 = i * a_dx;
    double x2 = (i + 1) * a_dx;
    double y1 = j * a_dy;
    double y2 = (j + 1) * a_dy;
    double divider = (x2 - x1) * (y2 - y1);
    double xx1 = x - x1;
    double xx2 = x2 - x;
    double yy1 = y - y1;
    double yy2 = y2 - y;
    double height = Q11 / divider * xx2 * yy2
                  + Q21 / divider * xx1 * yy2
                  + Q12 / divider * xx2 * yy1
                  + Q22 / divider * xx1 * yy1;
	return height;
}

/// reads level, waves, current data from file
void InteractiveWaves::setupArea(FileReader &fr, REAL dt, Initialization *pinit, REAL *ptime, std::string* pFile)
{
    m_Nx1 = fr.ReadInt();
    m_Nx2 = fr.ReadInt();
    m_dx1 = fr.ReadFloat();
    m_dx2 = fr.ReadFloat();
    m_h0_c = fr.ReadFloat();
    m_U1_c = fr.ReadFloat();
    m_U2_c = fr.ReadFloat();
    std::string spec = fr.ReadString();
    std::string waves = fr.ReadString();

    printf("%Gx%Gm = %d points\n", m_Nx1 * m_dx1, m_Nx2 * m_dx2, m_Nx1 * m_Nx2);

    // allocate the required arrays 0:Nx1+1,0:Nx2+1
    int sizeX = m_Nx1 + 3;
    int sizeY = m_Nx2 + 3;
    m_ph0 = new Array2D<REAL > (sizeX, sizeY);
    m_pU1 = new Array2D<REAL > (sizeX, sizeY);
    m_pU2 = new Array2D<REAL > (sizeX, sizeY);
    m_pwet = new Array2D<SHORT>(sizeX, sizeY);
    m_pzeta = new Array2D<REAL > (sizeX, sizeY);
    m_pzetam1 = new Array2D<REAL > (sizeX, sizeY);
    m_pphi = new Array2D<REAL > (sizeX, sizeY);
    m_pphim1 = new Array2D<REAL > (sizeX, sizeY);
    m_pzetaNet = new Array2D<REAL > (sizeX, sizeY);

    // fill the depth and mean velocity fields
    if (spec == std::string("-")) // use constant depth and mean velocities
    {
        *m_ph0 = m_h0_c;
        *m_pU1 = m_U1_c;
        *m_pU2 = m_U2_c;
        *m_pwet = 1;
        m_h0_min = (*m_ph0)[0][0];
        m_h0_max = (*m_ph0)[0][0];
    }
    else // read depth and velocities from file
    {
        FILE *fp = fopen(spec.c_str(), "r");
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
            {
                REAL a, b, c;
                bool d;
                ReadLine(fp, &a, &b, &c, &d);
                (*m_ph0)[i1][i2] = a;
                (*m_pU1)[i1][i2] = b;
                (*m_pU2)[i1][i2] = c;
                (*m_pwet)[i1][i2] = d ? 1 : 0;
                m_h0_min = std::min((*m_ph0)[i1][i2], m_h0_min);
                m_h0_max = std::max((*m_ph0)[i1][i2], m_h0_max);

                // maximize depth
                if (a > MAXIMUM_DEPTH)
                {
                    (*m_ph0)[i1][i2] = MAXIMUM_DEPTH;
                }
                else if (a < 0)
                {
                    (*m_ph0)[i1][i2] = 0;
                }
            }
        }
        fclose(fp);
        
#ifdef spiegelen
        // omschrijven in geval van spiegelbeeld
        fp = fopen("output.hu", "w");
        for (int i2 = m_Nx2 + 1; i2 >= 0; --i2)
        {
            for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
            {
                fprintf(fp, "%.6f %.6f %.6f %c\n",  (*m_ph0)[i1][i2], (*m_pU1)[i1][i2], 
                        (*m_pU2)[i1][i2], ((*m_pwet)[i1][i2] > 0 ? 'T' : 'F'));
            }
        }
        fclose(fp);
#endif
    }
    printf("max depth = %.2f min depth = %.2f\n", m_h0_max, m_h0_min);

    // fill the zeta and phi arrays
    if (waves == std::string("-") || waves == std::string("spectrum"))
        // initially everything is at rest or initialized via a wave spectrum
    {
        *pinit = waves == std::string("-") ? noInit : spectrumInit;
        *ptime = 0.0;
        *m_pzeta = 0.0;
        *m_pzetam1 = 0.0;
        *m_pphi = 0.0;
        *m_pphim1 = 0.0;
        *m_pzetaNet = 0.0; // nodig?
    }
    else // read waveheights and speed potentials from file
    {
        *pinit = fileInit;
        *pFile = waves;
    }

    // test of u1 en u2 verwaarloosd kunnen worden
    //    *m_pU1 = m_U1_c;
    //    *m_pU2 = m_U2_c;

    /* some numerical checks */
    
    // test if timestep is sufficiently small wrt grid size and depth
    double grid = std::min(m_dx1, m_dx2);
    REAL maxtimestep = grid / (2 * sqrt(m_g * m_h0_max));
    printf("h=%G dx=%G => dtmax=%G (dt=%G)\n", m_h0_max, grid, maxtimestep, dt);
    if (dt > maxtimestep)
    {
        fprintf(stderr, "Timestep too long (%f but must be < %f)\n", dt, maxtimestep);
        exit(-1);
    }
    
    REAL acc1 = 0.5 * m_g * FloatUtils::square(dt) / m_dx1;
    REAL acc2 = 0.5 * m_g * FloatUtils::square(dt) / m_dx2;
    REAL cfl1 = 0.0;
    REAL cfl2 = 0.0;
    REAL cfl1_max = 0.0;
    REAL cfl2_max = 0.0;
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            cfl1 = (fabs((*m_pU1)[i1][i2]) + sqrt(m_g * (*m_ph0)[i1][i2])) * dt / m_dx1;
            cfl2 = (fabs((*m_pU2)[i1][i2]) + sqrt(m_g * (*m_ph0)[i1][i2])) * dt / m_dx2;
            cfl1_max = std::max(cfl1, cfl1_max);
            cfl2_max = std::max(cfl2, cfl2_max);
        }
    }

    // check if timestep is not too large or spatial step too small
    // according to dt < dx / (2 * sqrt(g * h))
    if (FloatUtils::square(cfl1) + FloatUtils::square(cfl2) +
        FloatUtils::square(acc1) + FloatUtils::square(acc2) > 0.95)
    {
        printf("parameter problem\n");
        exit(-1);
    }

}


/// reads incoming waves specification from file
void InteractiveWaves::setupIncomingWaves(FileReader &fr)
{
    REAL pi = 4.0 * atan(1.0);
    WaveSpectrum waveType = static_cast<WaveSpectrum>(fr.ReadInt() + 1); // 0=deterministic, 1=jonswap in par files!
    m_wvin_spec_count = fr.ReadInt();
    int size = m_wvin_spec_count + 1;
    m_wvin_spec_ampl.resize(size);
    m_wvin_spec_phase.resize(size);
    m_wvin_spec_freq.resize(size);
    m_wvin_spec_dir.resize(size);
    if (waveType == deterministicWS)
    {
        for (int i = 1; i <= m_wvin_spec_count; ++i)
        {
            fr.ReadFloats(&m_wvin_spec_ampl[i], &m_wvin_spec_phase[i],
                          &m_wvin_spec_freq[i], &m_wvin_spec_dir[i]);
            // be sure the directions are in the range from 0 to 2*pi :
            m_wvin_spec_dir[i] = FloatUtils::modulo(m_wvin_spec_dir[i], 2.0 * pi);
            while (m_wvin_spec_dir[i] < 0)
            {
                m_wvin_spec_dir[i] += 2.0 * pi;
            }
        }
    }
    else if (waveType == jonswapWS)
    {
        REAL wvin_Hm0 = fr.ReadFloat();
        REAL wvin_gam_jon = fr.ReadFloat();
        REAL wvin_ome_m = fr.ReadFloat();
        REAL wvin_ome_max = fr.ReadFloat();
        REAL wvin_theta_m = fr.ReadFloat();
        REAL wvin_theta_std = fr.ReadFloat();
        REAL wvin_theta_max = fr.ReadFloat();

        /* compute random wave components */
        drawwaves(wvin_Hm0, wvin_gam_jon, wvin_ome_m, wvin_ome_max,
                  wvin_theta_m, wvin_theta_std, wvin_theta_max);
    }
    else
    {
        ASSERTFALSE(); // unknown incoming wave type in parameter file
    }
}


/// computes initial random waves
/** Draws random frequencies, phases and directions according to a Jonswap 
    spectrum and a truncated Gaussian directional distribution. */
void InteractiveWaves::drawwaves(REAL Hm0, ///< significant wave height [m]
                                 REAL gam_jon, ///< peak enhancement of Jonswap spectrum
                                 REAL ome_m, ///< mean angular frequency [rad/s]
                                 REAL ome_max, ///< maximum angular frequency [rad/s]
                                 REAL theta_m, ///< mean wave direction [rad]
                                 REAL theta_std, ///< standard deviation in the directional spreading function (0=unidir.) [rad]
                                 REAL theta_max)///< max. deviation from the mean direction [rad]
{
    REAL Sjon_max, Sjon_top, theta_dev, xa, pi, fac0, fac1, fac2, fac3, Sjon, sigma, ome_now;
    int j;
    bool not_found;

    pi = atan(1.0) * 4.0;

    Sjon_max = gam_jon * exp(-1.25);
    Sjon_top = 0.0;
    fac0 = sqrt(2.0 * (1.0 / 8.0 * FloatUtils::square(Hm0)) / m_wvin_spec_count);
    for (j = 1; j <= m_wvin_spec_count; ++j)
    {
        m_wvin_spec_ampl[j] = fac0;

        /* draw frequency from Jonswap spectrum */
        not_found = true;
        while (not_found)
        {
            /* draw a frequency */
            ome_now = ome_max * FloatUtils::randomNumber(false);

            /* compute Jonswap spectrum */
            xa = ome_now / ome_m;
            if (xa < 1.0)
            {
                sigma = 0.07;
            }
            else
            {
                sigma = 0.09;
            }
            fac1 = FloatUtils::power(xa, -5);
            fac2 = exp(-1.25 * FloatUtils::power(xa, -4));
            fac3 = FloatUtils::power(gam_jon, exp(-(FloatUtils::square(xa - 1)) / (2 * FloatUtils::square(sigma))));
            Sjon = fac1 * fac2 * fac3 / Sjon_max;
            Sjon_top = std::max(Sjon, Sjon_top);
            if (Sjon > 1.0)
            {
                ASSERTFALSE(); // wrong normalisation
            }
            if (FloatUtils::randomNumber(false) <= sqrt(Sjon))
            {
                not_found = false;
            }
        }
        m_wvin_spec_freq[j] = ome_now;

        /* draw phase */
        m_wvin_spec_phase[j] = 2 * pi * FloatUtils::randomNumber(false);

        /* draw wave propagation direction */
        if (theta_std > 0.0)
        {
            not_found = true;
            while (not_found)
            {
                theta_dev = sqrt(2.0) * theta_std * FloatUtils::randomNumber(true);
                if (fabs(theta_dev) <= theta_max)
                {
                    not_found = false;
                }
            }
        }
        else
        {
            theta_dev = 0.0;
        }
        m_wvin_spec_dir[j] = FloatUtils::modulo(theta_m + theta_dev, 2 * pi);
//        printf("%d %f %f %f %f\n", j - 1, m_wvin_spec_ampl[j], m_wvin_spec_freq[j], m_wvin_spec_phase[j], m_wvin_spec_dir[j]);
    }
} /* drawwaves */

// compute new net wave field (by adding pressure pulses to it) for use by own ships
void InteractiveWaves::correctForPulses(const Array2D<REAL>& pulseField)
{
    REAL invg = 1.0 / getGravity();
    ASSERT(pulseField.dim1() == m_pzetaNet->dim1(), "");
    ASSERT(pulseField.dim2() == m_pzetaNet->dim2(), "");

    int n1 = pulseField.dim1();
    int n2 = pulseField.dim2();
    for (int i = 0; i < n1; ++i)
    {
        for (int j = 0; j < n2; ++j)
	{
            (*m_pzetaNet)[i][j] = (*m_pzeta)[i][j] + pulseField[i][j] * invg;
#ifdef PRINTDEBUG
            if (fabs((*m_pzetaNet)[i][j]) > 100)
            {
                printf("%G = %G + %G * %G\n", (*m_pzetaNet)[i][j], (*m_pzeta)[i][j], pulseField[i][j], invg);
            }
#endif
        }
    }
}
